#### Figure 2

# Clear
rm(list=ls())

# Necessary packages
library(dplyr)
library(haven)

# Working directory
setwd("~/Dropbox/Paraguay/1_data/7_replication_dataverse")

# Load land data
tmh_final <- read_dta("data/land.dta")


# Collapse to year level
tmh_yr <- tmh_final %>%
  group_by(year) %>%
  summarise(
    tmh.ha = sum(tmh_ha),          # Sum total hectares
    tmh.n = n_distinct(ID)         # Count unique families 
  ) %>%
  ungroup()

# Convert hectares to thousands and round
tmh_yr$tmh.ha <- tmh_yr$tmh.ha / 1000
tmh_yr$tmh.ha <- round(tmh_yr$tmh.ha, 1)

# figure2a: Hectares time series
pdf(file = "figures/figure2a.pdf")
plot(tmh_yr$tmh.ha ~ tmh_yr$year, 
     type = "h", 
     lty = 1,
     lwd = 8.5, 
     cex.lab = 1.15,
     col = "darkgrey",
     axes = FALSE,
     xlim = c(1950, 2005),
     ylim = c(0, 900),
     xlab = "Year", 
     ylab = "Hectares (in thousands)", 
     frame.plot = FALSE)
abline(v = 1989, col = "darkred", lty = 2, lwd = 2)
axis(side = 1, at = seq(1950, 2005, 5), cex.axis = 1.15)
axis(side = 2, at = seq(0, 900, 300), cex.axis = 1.15)
dev.off()

# figure2b: Families time series
pdf(file = "figures/figure2b.pdf")
plot(tmh_yr$tmh.n ~ tmh_yr$year, 
     type = "h", 
     lty = 1,
     lwd = 8.5, 
     cex.lab = 1.15,
     col = "darkgrey",
     axes = FALSE,
     xlim = c(1950, 2005),
     ylim = c(0, 200),
     xlab = "Year", 
     ylab = "Families", 
     frame.plot = FALSE)
abline(v = 1989, col = "darkred", lty = 2, lwd = 2)
axis(side = 1, at = seq(1950, 2005, 5), cex.axis = 1.15)
axis(side = 2, at = seq(0, 200, 50), cex.axis = 1.15)
dev.off()
